CO2 Na França
Grupo:
No presente projeto, iremos, a partir de base de dados disponibilizada em aula, gerar um modelo linear para a previsão de CO2 em KT. Os resultados obtidos foram positivos e seguiram todos os requisitos necessários para a validade da regressão:
O modelo foi desenvolvido, no âmbito de nosso grupo, para os dados relativos à França. Nas próximas seções, iremos descrever cada um dos passos realizados, desde a análise exploratória de dados até os resultados finais, em formato de storytelling.
Iremos iniciar o trabalho verificando os dados que temos para a realização da análise. Para isso, seguiremos a seguinte sequência de passos:
setwd('~/Área de Trabalho/stu/MBA - FGV/Estatística Avançada/Trabalhos/Trabalho-Final')
df_in <- read_excel('./FRA_Country_en_excel_v2.xls')
## New names:
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * `` -> ...6
## * `` -> ...7
## * … and 53 more problems
df_in[, 1:3] %>% head
Os dados parecem extremamente desorganizados. Precisamos, primeiramente, eliminar as duas primeiras linhas e tornar a terceira linha o header do Data Frame:
colnames(df_in) <- df_in[3,] %>% unlist
df_in <- df_in[-c(1, 2, 3),]
df_in[, 1:3] %>% head
Sabemos que o país é a França, então podemos eliminar o nome e o código do país. O nome do indicador será utilizado para se referir à variável e o nome do indicador será eliminado da tabela. Um segundo Data Frame relacionando os códigos dos indicadores aos nomes será criado e funcionará como uma tabela de metadados:
df_in_meta <- df_in %>% select(`Indicator Name`, `Indicator Code`)
df_in_meta %>% head
Eliminando os metadados e as informações redundantes de país da tabela de dados:
df_in <- df_in %>% select(-`Country Name`, -`Country Code`, -`Indicator Name`)
df_in[, 1:3] %>% head
Podemos notar que existem muitos indicadores com um altíssimo número de dados perdidos. Iremos então verificar em uma tabela a quantidade de dados perdidos por indicador, para descobrir quais indicadores não serão levados em consideração na análise:
df_in$`N. Missing` <- apply(df_in, 1, function(x) (is.na(x) %>% sum))
df_in$`Perc. Missing` <- df_in$`N. Missing` / (ncol(df_in) - 2)
df_in_missing <- df_in %>%
select(`Indicator Code`, `N. Missing`, `Perc. Missing`) %>%
arrange(desc(`N. Missing`))
Escrevendo a a saída em uma tabela Shiny formatada e em um gráfico de barra com o percentual de dados perdidos:
ggplot(data=df_in_missing, aes(x=`N. Missing`)) +
geom_histogram(color='darkgreen', fill='white') + theme_minimal() +
labs(x='Quantidade de dados perdidos',
y='Freq. Absoluta',
title='Missing Data',
subtitle='Analisando quantidade de dados perdidos por tipo de variável')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
df_in_missing
Se tomarmos um limiar igual a 5, admitiremos um valor de perda máximo de dados igual a a aproximadamente 10% e teremos no mínimo 50 pontos de amostra para gerar o modelo:
ggplot(data=df_in_missing %>% filter(`N. Missing` <= 5), aes(x=`N. Missing`)) +
geom_histogram(color='darkgreen', fill='white') + theme_minimal() +
labs(x='Quantidade de dados perdidos',
y='Freq. Absoluta',
title='Missing Data',
subtitle='Filtrando gráfico anterior para regiões com menos de 6 dados perdidos')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
df_in_missing %>% filter(`N. Missing` <= 5)
De acordo com o histograma, ainda que tal critério rigoroso seja adotado, uma grande quantidade de indicadores continua disponível: reduzimos o escopo da análise de 1200 indicadores para 200 indicadores e, ainda assim, temos um grande volume de dados.
Podemos, então, plotar conjuntamente os diagramas para verificar, de maneira superficial, o aspecto das curvas que serão submetidas a análise. Para facilitar tal tarefa, iremos transpor a tabela filtrada:
df_in <- df_in %>% filter(`N. Missing` <= 5)
data_matrix <- df_in[1:nrow(df_in), 2:(ncol(df_in) - 2)]
df_in_t <- t(data_matrix) %>% as.data.frame
colnames(df_in_t) <- df_in$`Indicator Code`
df_in_t$Year <- colnames(df_in)[2:(ncol(df_in) - 2)]
df_in_t[,1:3] %>% head
Iremos eliminar o ano de \(2015\) da análise pois foi observado que a coluna de tal ano não possui nenhum valor para nenhuma variável.
Plotando a evolução de cada uma das variáveis em função do tempo:
df_in_p <- df_in_t %>% gather('Indicator', 'Value', -Year)
## Warning: attributes are not identical across measure variables;
## they will be dropped
df_in_p <- df_in_p %>% filter(Year != '2015')
df_in_p <- df_in_p %>% transform(Year=as.numeric(Year)) %>% transform(Value=as.numeric(Value))
ggplot(df_in_p, aes(x=Year, y=Value, color=Indicator)) + geom_line() +
theme(legend.position='none', panel.background = element_blank(),
panel.grid.major = element_line(colour = 'gray'),
panel.grid.minor = element_line(colour = 'gray')) +
labs(x='Ano',
y='Valor',
title='Evolução das Variáveis Escolhidas',
subtitle='Uma primeira análise das séries temporais')
## Warning: Removed 148 rows containing missing values (geom_path).
Podemos observar que existem dados perdidos devidos aos “Warnings” recebidos. Na próxima seção iremos corrigir esse problema.
A função “approx” da biblioteca zoo permite que realizemos interpolação linear para preencher os dados faltantes de cada uma das séries. Esse procedimento é mais coerente que simplesmente copiar os dados do passado ou do futuro para completar as células iguais a NA.
Esse procedimento não considera os pontos extremos com valores iguais a NA, razão pela qual, após a interpolação, escrevemos uma rotina para preencher os valores extremos com o próximo valor nao nulo ou o valor não nulo imediatamente anterior.
indicators_list <- colnames(df_in_t %>% select(-Year))
for (curr_indicator in indicators_list) {
curr_list <- df_in_t[curr_indicator]
curr_list <- na.approx(curr_list, na.rm = FALSE)
non_NA_index <- which(!is.na(curr_list))
first_non_NA <- min(non_NA_index)
last_non_NA <- max(non_NA_index)
if (first_non_NA > 1){
curr_list[1:(first_non_NA - 1)] = curr_list[first_non_NA]
}
if (last_non_NA < length(curr_list)) {
curr_list[(last_non_NA + 1):length(curr_list)] = curr_list[last_non_NA]
}
df_in_t[curr_indicator] = curr_list
}
Finalmente, podemos plotar novamente os gráficos:
df_in_p <- df_in_t %>% gather('Indicator', 'Value', -Year)
## Warning: attributes are not identical across measure variables;
## they will be dropped
df_in_p[is.na(df_in_p)] <- 0
df_in_p <- df_in_p %>% transform(Year=as.numeric(Year)) %>% transform(Value=as.numeric(Value))
ggplotly(ggplot(df_in_p, aes(x=Year, y=Value, color=Indicator)) + geom_line() +
theme(legend.position='none', panel.background = element_blank(),
panel.grid.major = element_line(colour = 'gray'),
panel.grid.minor = element_line(colour = 'gray')) +
labs(x='Ano (Variáveis reativas: encostar mouse para verificar)',
y='Valor',
title='Variáveis Corrigidas'))
Não há alterações visíveis pois os dados perdidos se concentraram majoritariamente em posições extremas das seŕies. Em todo caso, os warnings de dados faltantes não se encontram mais presentes e podemos prosseguir para a próxima etapa: analisar a correlação e o relacionamento entre as variáveis.
Podemos plotar a correlação entre as variáveis em um mapa de calor. É inviável escrever o nome de cada uma das variáveis na matriz, entretanto, podemos, ao menos, observar o aspecto de tal mapa.
cor_mat <- df_in_t %>% select(-Year) %>% cor
ggcorrplot(cor_mat, tl.cex=0) +
labs(title='Matriz de correlação como Mapa de Calor',
subtitle='Identificando se há multicolinearidade')
Aparentemente, os grupos de variáveis são fortemente relacionados entre si. Isso indica que podemos agrupar as variáveis em um número muito mais compacto de categorias sem que isso signifique prejuízo em nossa análise.
Vejamos a correlação entre cada uma das variáveis com o dado em estudo na próxima subseção.
Antes de agruparmos as variáveis em categorias correlacionadas, vamos dar um “Zoom” na análise da seção anterior e nos ater à variável em estudo. Trata-se da emissão de CO² em KT, representada pela variável EN.ATM.CO2E.KT.
Podemos visualizar as correlações de cada uma das demais variáveis por meio de um gráfico de barra.
x_plot <- rownames(cor_mat)
y_plot <- cor_mat[, 'EN.ATM.CO2E.KT']
df_plot <- data.frame(X=x_plot, Y=y_plot)
ggplotly(ggplot(data=df_plot, aes(x=reorder(x_plot, -abs(y_plot)), y=abs(y_plot))) +
geom_bar(stat='identity', color='darkgreen', fill='white') +
theme_minimal() + theme(axis.text.x=element_blank()) +
labs(x='Nome da Variável (Gráfico reativo: encostar o mouse para verificar)',
y='Correlação com Emissão de CO2 (KT)',
title='Correlações com Variável Alvo'))
Plotando o valor absoluto da correlação:
ggplotly(ggplot(data=df_plot, aes(x=reorder(x_plot, -abs(y_plot)), y=abs(y_plot))) +
geom_bar(stat='identity', color='darkgreen', fill='white') +
theme_minimal() + theme(axis.text.x=element_blank()) +
labs(x='Nome da Variável (Gráfico reativo: encostar o mouse para verificar)',
y='Correlação com Emissão de CO2 (KT)',
title='Correlações com Variável Alvo'))
Não é, a princípio, necessário saber o nome de cada uma das variáveis representadas por cada uma das barras. O que sabemos, por hora, é que a correlação varia praticamente de forma contínua e todo tipo de correlação ocorre no grande conjunto de dados que temos disponíveis.
Assim, na próxima seção, iremos tentar agrupar as variáveis em grupos correlacionados no intuito de:
Antes de efetuarmos o agrupamento, vamos retirar de nossa tabela a variável que queremos explicar>
df_out <- df_in_t[, 'EN.ATM.CO2E.KT']
df_in_t <- df_in_t %>% select(-EN.ATM.CO2E.KT)
As técnicas de aprendizado não supervisionado consistem na categorização de variáveis por meio de algoritmos de Cluster. Esse tipo de procedimento pode ser aplicado no agrupamento de variáveis fortemente correlacionadas.
Para isso, podemos imaginar que quanto mais correlacionadas duas variáveis forem, mais próximas elas estarão no hiperespaço de variáveis explicativas. Ou seja: podemos adotar uma métrica \(\mathcal{F}(Cor(X, Y)) = \mathcal{D}(X, Y)\) que corresponde à distância entre os vetores \(X\) e \(Y\).
Essa métrica é tal que, quanto mais correlacionadas forem as variáveis, menor será a distância entre elas, de tal forma que variáveis \(100\%\) correlacionadas serão separadas por uma distância igual a zero.
Além disso, a correlação será avaliada em valor absoluto neste ponto. Isso ocorre porque variáveis com correlações próximas a \(-1\) também podem ser consideradas extremamente próximas, ainda que as variações de cada uma delas ocorram com sinais trocados.
Definiremos a função em questão como: \(\mathcal{F}(Cor(X, Y)) = \mathcal{D}(X, Y) = 1 - |Cor(X, Y)|\).
A correlação entre cada par de variáveis formará uma matriz de ordem \(N \times N\), onde \(N\) é a quantidade de variáveis presentes na base de dados (em torno de 200 variáveis).
Assim, essa matriz terá a forma: \(\mathcal{D}_K(X_i, X_j) = 1 - |Cor(X_i, X_j)| = (d)_{ij} = D_{N \times N}\)
Os pontos serão agrupados por meio da técnica de clusterização hierárquica, considerando-se a matriz de distância definida acima: \(\mathcal{F_K}(Cor(X, Y)) = \mathcal{D}(X, Y) = 1 - |Cor(X, Y)|\).
A clusterização hierárquica não pede que o usuário determine o número de clusters desejado. Ao contrário, o usuário pode visualizar um dendograma com diferentes níveis de agrupamento e escolher aquele que melhor se aplica à necessidade do problema:
hclust_obj <- (1 - abs(cor_mat)) %>% as.dist %>% hclust
dend <- hclust_obj %>% as.dendrogram
dend %>% dendextend::set('labels_color', 'white') %>% plot(xaxt='n')
Precisamos escolher uma altura na qual a árvore será cortada. Para facilitar essa análise, plotemos o número de clusters em função da altura de corte:
height <- seq(from=0, to=1, by=0.001)
n_clusters <- sapply(height, function(X)(cutree(dend, h=X) %>% unique %>% length))
df_hclust <- data.frame(Height=height, N.Clusters=n_clusters)
ggplotly(ggplot(df_hclust, aes(x=Height, y=N.Clusters)) + theme_minimal() +
geom_line(stat='identity', color='darkgreen'))
O gráfico plotado é reativo, podemos encostar o mouse sobre ele e verificar, manualmente, os valores de \(X\) (altura) e \(Y\) (número de clusters) para cada uma das barras. Podemos observar a existência de um platô quando a altura do algoritmo é igual a \(0.2\). Nesse caso, deixamos de trabalhar com \(200\) variáveis e passamos a operar com \(40\) grupos diferentes.
Iremos então acessar o nosso dicionário de metadados e atribuir a cada uma das variáveis o respectivo identificador do cluster.
opt_cut <- cutree(dend, h=0.2)
df_in_clusters <- data.frame('Indicator.Code'=names(opt_cut),
'Cluster.Index'=as.vector(opt_cut))
df_in_clusters <- df_in_clusters %>%
inner_join(df_in_meta, by=(c('Indicator.Code'='Indicator Code')))
## Warning: Column `Indicator.Code`/`Indicator Code` joining factor and
## character vector, coercing into character vector
df_in_clusters %>% head
write.csv(df_in_clusters, './Agrupamento_Variaveis.csv')
Finalmente, podemos observar o arquivo de saída e: * Descrever cada um dos grupos * Interpretar e caracterizar cada cluster * Identificar potenciais correlações espúrias
Isso será feito na próxima seção.
Uma tabela caracterizando cada um dos clusteres foi criada. Ela será impressa nesta seção para consulta e comparação com as análises que serão feitas daqui em diante. As descrições foram feitas relacionando-se cada grupo de variáveis dentro de um mesmo cluster e interpretando o que cada agrupamento representa.
df_grupos <- read_excel('./GRUPOS.xlsx')
df_grupos
Podemos avaliar a qualidade dos clusters verificando a distribuição das correlações para agrupamento. Queremos grupos com distribuições de correlações o mais próximas de \(100\%\) que for possível.
get_corr_distrib <- function(cluster_index) {
cluster_table <- df_in_clusters %>% filter(Cluster.Index == cluster_index)
vars_list <- cluster_table[['Indicator.Code']]
n_cluster_vars <- length(vars_list)
if (n_cluster_vars == 1) {
return(NA)
}
corr_vec <- c()
for (i in 2:n_cluster_vars) {
for (j in 1:(i - 1)) {
corr_vec <- c(corr_vec, cor_mat[vars_list[[i]], vars_list[[j]]] %>% abs)
}
}
return(corr_vec)
}
corr_vec <- c()
cluster_vec <- c()
for (k in (df_in_clusters[['Cluster.Index']] %>% unique)) {
new_rows <- get_corr_distrib(k)
corr_vec <- c(corr_vec, new_rows)
cluster_vec <- c(cluster_vec, rep(k %>% as.character, times=new_rows %>% length))
}
df_corr_per_cluster <- data.frame(Correlation=corr_vec, Cluster.Index=cluster_vec) %>%
mutate(Cluster.Index=as.numeric(Cluster.Index)) %>%
inner_join(df_grupos, by=c('Cluster.Index'='Cluster.Index'))
ggplot(df_corr_per_cluster) +
geom_density(aes(x=Correlation, fill=Cluster.Cod), alpha=.1) +
scale_x_continuous(limits=c(.7, 1)) + theme_minimal() +
labs(x='Correlação', y='Densidade por Cluster',
title='Distribuição da correlação',
subtitle='Por pares de variáveis dentro de um mesmo grupo / por grupo')
## Warning: Removed 17 rows containing non-finite values (stat_density).
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
Podemos ainda plotar a distribuição global, sem segregação por cluster:
ggplot(df_corr_per_cluster, aes(x=Correlation)) + geom_density() +
scale_x_continuous(limits=c(.7, 1)) + theme_minimal() +
labs(x='Correlação (Variáveis em um mesmo cluster)', y='Densidade',
title='Distribuição da correlação',
subtitle='Por pares de variáveis dentro de um mesmo grupo / geral')
## Warning: Removed 17 rows containing non-finite values (stat_density).
Podemos então observar que, além de reduzir o número de variáveis de \(200\) para \(32\), garantimos que dentro de um cluster teremos correlações iguais a, no mínimo, \(70\%\). Assim, podemos prosseguir nossa análise, focando apenas em variáveis agregadas por grupos. A agregação será realizada em duas etapas:
Variáveis com correlação maior ou igual a \(99\%\) podem ser consideradas como sendo a mesma variável. Assim, iremos em uma primeira etapa retirar a média de grupos de variáveis com tal correlação fortíssima.
Em seguida, agruparemos o cluster realizando, também, uma média aritimética. Nesta etapa, os erros gerados por cada variável também serão submetidos a tal média e isso contribui na normalização do erro (pela lei do limite central), na redução da variância do erro (pois o desvio padrão se reduz em um fator igual a \(\sqrt{N_{Cluster}}\)) e na eliminação do problema da multicolinearidade sem eliminarmos qualquer variável da análise.
O agrupamento por média deverá ser efetuado sobre as variáveis normalizadas, para que o efeito da média não seja comparado a uma média ponderada na qual variáveis com unidades de maiores ordens de grandeza prevalecem.
df_in_p_sd <- df_in_p
df_in_p_sd <- df_in_p_sd %>%
group_by(Indicator) %>%
summarise_at(vars(Value), sd) %>%
rename(Sd = Value) %>%
as.data.frame
df_in_p_mean <- df_in_p
df_in_p_mean <- df_in_p_mean %>%
group_by(Indicator) %>%
summarise_at(vars(Value), mean) %>%
rename(Mean = Value) %>%
as.data.frame
df_in_p_grupos <- df_in_p %>%
mutate(Indicator = as.character(Indicator)) %>%
inner_join(df_in_p_sd, by=c('Indicator'='Indicator')) %>%
inner_join(df_in_p_mean, by=c('Indicator' = 'Indicator')) %>%
mutate(Value = (Value - Mean) / Sd) %>%
inner_join(df_in_clusters, by=c('Indicator' = 'Indicator.Code')) %>%
mutate(Cluster.Index = as.numeric(Cluster.Index)) %>%
inner_join(df_grupos, by=c('Cluster.Index' = 'Cluster.Index')) %>%
select('Cluster.Index', 'Cluster.Cod', 'Cluster.Description', 'Year', 'Value') %>%
group_by(Cluster.Cod, Year) %>%
summarise_at(vars(Value), mean) %>%
as.data.frame
df_in_p_grupos %>% head
Há correlação entre grupos de um mesmo cluster? É necessário fazer essa verificação final para conferir se o problema da multicolinearidade pode ocorrer ainda. Iniciemos plotando a evolução do valor médio encontrado em cada um dos clusters:
clusters_list <- df_in_p_grupos[['Cluster.Cod']] %>% unique
ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[1:6]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
subtitle='Variáveis 1 a 6')
ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[7:12]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
subtitle='Variáveis 7 a 12')
ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[13:18]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
subtitle='Variáveis 13 a 18')
ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[19:24]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
subtitle='Variáveis 19 a 24')
ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[25:30]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
subtitle='Variáveis 25 a 30')
ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[31:36]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
subtitle='Variáveis 31 a 36')
ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[37:39]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
subtitle='Variáveis 37 a 39')
Para nos certificarmos de maneira exata, verifiquemos o correlograma. Para isso, precisaremos retirar a tabela utilizada nos gráficos anteriores da forma de pivô.
list_years <- df_in_p_grupos[['Year']] %>% unique
list_clusters <- df_in_p_grupos[['Cluster.Cod']] %>% unique
nrows_mat <- length(list_years)
ncols_mat <- length(list_clusters)
df_in_t_grupos <- zeros(nrows_mat, ncols_mat)
rownames(df_in_t_grupos) <- list_years
colnames(df_in_t_grupos) <- list_clusters
for (curr_year in list_years) {
for (curr_cluster in list_clusters) {
df_in_t_grupos[[curr_year %>% as.character, curr_cluster]] <-
(df_in_p_grupos %>%
filter(Year == curr_year &
Cluster.Cod == curr_cluster))$Value[[1]]
}
}
df_in_t_grupos <- df_in_t_grupos %>% as.data.frame
mat_cor_grupos <- cor(df_in_t_grupos) %>% abs
mask <- zeros(ncols_mat) > 0
for (i in 1:ncols_mat) {
mat_cor_grupos[[i, i]] <- 0
mask[[i]] <- (any(abs(mat_cor_grupos[i,]) > 0.78))
mat_cor_grupos[[i, i]] <- 1
}
mat_cor_grupos[mask, mask] %>% corrplot
Ainda temos um grande problema de multicolinearidade entre os clusters. Para resolvermos isso, podemos recorrer à técnica da análise de componentes principais.
Poderíamos, desde o início, ter realizado essa técnica, sem o intermédio de qualquer algoritmo de clusterização. No entanto, analisar a influência de um número reduzido de clusters em cada um dos componentes principais é mais fácil e garante uma maior interpretabilidade ao modelo.
pr_comp <- prcomp(df_in_t_grupos)
list_var <- (pr_comp$sdev ** 2) / (sum(pr_comp$sdev ** 2))
list_cum <- cumsum(list_var)
df_pr <- data.frame(Perc.Var=list_var,
Cum.Perc=list_cum,
PC.Index=1:length(list_var))
ggplot(df_pr, aes(x=PC.Index)) +
geom_bar(aes(y=Perc.Var), stat='identity', color='darkgreen', fill='white') +
geom_line(aes(y=Cum.Perc)) +
scale_x_continuous(breaks = seq(0, 40, len=10)) +
theme_minimal() +
labs(x='Índice da Componente Principal',
y='Percentual de Importância da Componente',
title='Gráfico de Pareto: Componentes Principais',
subtitle='Importância relativa (barras) e importância cumulativa (linha)')
pr_comp %>% summary
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.1518 2.1739 2.005 1.46340 1.22961 1.10821 1.0430
## Proportion of Variance 0.3361 0.1599 0.136 0.07245 0.05115 0.04155 0.0368
## Cumulative Proportion 0.3361 0.4960 0.632 0.70441 0.75556 0.79711 0.8339
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.9650 0.87051 0.79517 0.68265 0.58170 0.56293
## Proportion of Variance 0.0315 0.02564 0.02139 0.01577 0.01145 0.01072
## Cumulative Proportion 0.8654 0.89105 0.91244 0.92821 0.93965 0.95037
## PC14 PC15 PC16 PC17 PC18 PC19
## Standard deviation 0.49146 0.4582 0.4384 0.40714 0.38020 0.32840
## Proportion of Variance 0.00817 0.0071 0.0065 0.00561 0.00489 0.00365
## Cumulative Proportion 0.95855 0.9657 0.9721 0.97776 0.98265 0.98630
## PC20 PC21 PC22 PC23 PC24 PC25
## Standard deviation 0.28852 0.26693 0.24880 0.20955 0.19867 0.16624
## Proportion of Variance 0.00282 0.00241 0.00209 0.00149 0.00134 0.00093
## Cumulative Proportion 0.98911 0.99152 0.99362 0.99510 0.99644 0.99737
## PC26 PC27 PC28 PC29 PC30 PC31
## Standard deviation 0.14316 0.11484 0.10626 0.09295 0.08766 0.06892
## Proportion of Variance 0.00069 0.00045 0.00038 0.00029 0.00026 0.00016
## Cumulative Proportion 0.99807 0.99851 0.99890 0.99919 0.99945 0.99961
## PC32 PC33 PC34 PC35 PC36 PC37
## Standard deviation 0.06533 0.05669 0.04332 0.03426 0.02667 0.01283
## Proportion of Variance 0.00014 0.00011 0.00006 0.00004 0.00002 0.00001
## Cumulative Proportion 0.99975 0.99986 0.99992 0.99996 0.99999 0.99999
## PC38 PC39
## Standard deviation 0.01123 0.007274
## Proportion of Variance 0.00000 0.000000
## Cumulative Proportion 1.00000 1.000000
Podemos observar que grande parte das variações nas séries pode ser explicada por um pequeno número de componentes principais (conseguimos explicar 95% das variações a partir de 12 componentes). Assim, esperamos que, na etapa de modelagem, consigamos criar um modelo com poucas variáveis explicativas e buscaremos eliminá-las a partir de duas etapas:
Em cada instante de tempo, o vetor \(\mathcal{C} = (C_1, C_2,..., C_{39})\) formado pelos clusteres será multiplicado pela matriz de rotação encontrada na análise de componentes principais e, assim, encontraremos novas séries temporais sobre as quais serão efetuadas as análises de regressão.
Aplicando a transformação em nossa matriz de clusters, obtemos o termo “x” da saída da função prcomp.
df_pca <- pr_comp$x %>% as.data.frame
df_pca$Y <- df_out
Finalmente, podemos iniciar a modelagem do regressor linear.
Ajustando o modelo, nas componentes principais:
mod = lm(Y~., data=df_pca)
summary(mod)
##
## Call:
## lm(formula = Y ~ ., data = df_pca)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6328.6 -1902.9 -338.7 1546.3 6104.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 391149.2 674.6 579.861 < 2e-16 ***
## PC1 7948.6 216.0 36.807 < 2e-16 ***
## PC2 1689.9 313.1 5.397 5.93e-05 ***
## PC3 -20651.3 339.5 -60.832 < 2e-16 ***
## PC4 -19536.4 465.1 -42.003 < 2e-16 ***
## PC5 -5201.2 553.6 -9.396 6.49e-08 ***
## PC6 8640.3 614.2 14.068 1.99e-10 ***
## PC7 419.2 652.6 0.642 0.529784
## PC8 733.6 705.3 1.040 0.313762
## PC9 -2792.0 781.9 -3.571 0.002552 **
## PC10 -889.0 856.0 -1.039 0.314453
## PC11 -7610.6 997.1 -7.633 1.01e-06 ***
## PC12 -4530.4 1170.1 -3.872 0.001352 **
## PC13 8178.5 1209.2 6.764 4.55e-06 ***
## PC14 -6912.5 1385.0 -4.991 0.000133 ***
## PC15 6402.1 1485.6 4.310 0.000540 ***
## PC16 2178.5 1552.6 1.403 0.179703
## PC17 -3261.5 1671.8 -1.951 0.068809 .
## PC18 -3827.9 1790.3 -2.138 0.048277 *
## PC19 1670.9 2072.6 0.806 0.431963
## PC20 -3589.7 2359.1 -1.522 0.147613
## PC21 15742.9 2550.0 6.174 1.34e-05 ***
## PC22 9454.3 2735.8 3.456 0.003254 **
## PC23 15197.6 3248.2 4.679 0.000252 ***
## PC24 16195.3 3426.1 4.727 0.000228 ***
## PC25 6123.6 4094.4 1.496 0.154217
## PC26 13572.9 4754.5 2.855 0.011468 *
## PC27 1439.1 5927.0 0.243 0.811247
## PC28 -6166.5 6405.3 -0.963 0.350020
## PC29 41051.8 7323.0 5.606 3.94e-05 ***
## PC30 12587.3 7764.7 1.621 0.124535
## PC31 37646.4 9876.8 3.812 0.001535 **
## PC32 -9955.8 10418.3 -0.956 0.353486
## PC33 -32279.3 12006.6 -2.688 0.016149 *
## PC34 15625.8 15713.2 0.994 0.334811
## PC35 20042.1 19866.4 1.009 0.328062
## PC36 14107.2 25520.1 0.553 0.588049
## PC37 102290.8 53069.9 1.927 0.071866 .
## PC38 -108663.4 60589.4 -1.793 0.091824 .
## PC39 -78836.9 93575.7 -0.842 0.411929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5048 on 16 degrees of freedom
## Multiple R-squared: 0.9979, Adjusted R-squared: 0.9927
## F-statistic: 192.3 on 39 and 16 DF, p-value: 7.613e-16
Observamos que para diversas variáveis não podemos rejeitar a hipótese \(\mathcal{H}_0\) de nulidade do coeficiente. Assim, iremos tentar resolver o problema eliminando as variáveis em sequência adequada a aumentar a verossimilhança do modelo por meio da função “step”, utilizando-se o parâmetro “direction = backward”:
mod_ajustado <- step(mod, direction = 'backward')
## Start: AIC=964.84
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 +
## PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 +
## PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC27 + PC28 +
## PC29 + PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC36 + PC37 +
## PC38 + PC39
##
## Df Sum of Sq RSS AIC
## - PC27 1 1.5022e+06 4.0921e+08 963.05
## - PC36 1 7.7865e+06 4.1549e+08 963.90
## - PC7 1 1.0512e+07 4.1822e+08 964.26
## <none> 4.0770e+08 964.84
## - PC19 1 1.6561e+07 4.2427e+08 965.07
## - PC39 1 1.8087e+07 4.2579e+08 965.27
## - PC32 1 2.3269e+07 4.3097e+08 965.95
## - PC28 1 2.3616e+07 4.3132e+08 965.99
## - PC34 1 2.5199e+07 4.3290e+08 966.20
## - PC35 1 2.5934e+07 4.3364e+08 966.29
## - PC10 1 2.7484e+07 4.3519e+08 966.49
## - PC8 1 2.7565e+07 4.3527e+08 966.50
## - PC16 1 5.0163e+07 4.5787e+08 969.34
## - PC25 1 5.6999e+07 4.6470e+08 970.17
## - PC20 1 5.9001e+07 4.6671e+08 970.41
## - PC30 1 6.6964e+07 4.7467e+08 971.36
## - PC38 1 8.1959e+07 4.8966e+08 973.10
## - PC37 1 9.4668e+07 5.0237e+08 974.53
## - PC17 1 9.6981e+07 5.0469e+08 974.79
## - PC18 1 1.1650e+08 5.2420e+08 976.91
## - PC33 1 1.8418e+08 5.9188e+08 983.71
## - PC26 1 2.0766e+08 6.1536e+08 985.89
## - PC22 1 3.0431e+08 7.1201e+08 994.06
## - PC9 1 3.2490e+08 7.3261e+08 995.66
## - PC31 1 3.7020e+08 7.7790e+08 999.02
## - PC12 1 3.8197e+08 7.8968e+08 999.86
## - PC15 1 4.7324e+08 8.8094e+08 1005.98
## - PC23 1 5.5779e+08 9.6550e+08 1011.12
## - PC24 1 5.6938e+08 9.7708e+08 1011.78
## - PC14 1 6.3476e+08 1.0425e+09 1015.41
## - PC2 1 7.4226e+08 1.1500e+09 1020.91
## - PC29 1 8.0078e+08 1.2085e+09 1023.69
## - PC21 1 9.7120e+08 1.3789e+09 1031.08
## - PC13 1 1.1658e+09 1.5735e+09 1038.47
## - PC11 1 1.4845e+09 1.8922e+09 1048.80
## - PC5 1 2.2496e+09 2.6573e+09 1067.81
## - PC6 1 5.0428e+09 5.4505e+09 1108.04
## - PC1 1 3.4521e+10 3.4928e+10 1212.07
## - PC4 1 4.4955e+10 4.5363e+10 1226.71
## - PC3 1 9.4295e+10 9.4702e+10 1267.92
##
## Step: AIC=963.05
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 +
## PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 +
## PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + PC29 +
## PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC36 + PC37 + PC38 +
## PC39
##
## Df Sum of Sq RSS AIC
## - PC36 1 7.7865e+06 4.1699e+08 962.10
## - PC7 1 1.0512e+07 4.1972e+08 962.47
## <none> 4.0921e+08 963.05
## - PC19 1 1.6561e+07 4.2577e+08 963.27
## - PC39 1 1.8087e+07 4.2729e+08 963.47
## - PC32 1 2.3269e+07 4.3248e+08 964.14
## - PC28 1 2.3616e+07 4.3282e+08 964.19
## - PC34 1 2.5199e+07 4.3441e+08 964.39
## - PC35 1 2.5934e+07 4.3514e+08 964.49
## - PC10 1 2.7484e+07 4.3669e+08 964.69
## - PC8 1 2.7565e+07 4.3677e+08 964.70
## - PC16 1 5.0163e+07 4.5937e+08 967.52
## - PC25 1 5.6999e+07 4.6621e+08 968.35
## - PC20 1 5.9001e+07 4.6821e+08 968.59
## - PC30 1 6.6964e+07 4.7617e+08 969.53
## - PC38 1 8.1959e+07 4.9117e+08 971.27
## - PC37 1 9.4668e+07 5.0387e+08 972.70
## - PC17 1 9.6981e+07 5.0619e+08 972.96
## - PC18 1 1.1650e+08 5.2570e+08 975.07
## - PC33 1 1.8418e+08 5.9338e+08 981.86
## - PC26 1 2.0766e+08 6.1687e+08 984.03
## - PC22 1 3.0431e+08 7.1351e+08 992.18
## - PC9 1 3.2490e+08 7.3411e+08 993.77
## - PC31 1 3.7020e+08 7.7941e+08 997.13
## - PC12 1 3.8197e+08 7.9118e+08 997.97
## - PC15 1 4.7324e+08 8.8245e+08 1004.08
## - PC23 1 5.5779e+08 9.6700e+08 1009.20
## - PC24 1 5.6938e+08 9.7858e+08 1009.87
## - PC14 1 6.3476e+08 1.0440e+09 1013.49
## - PC2 1 7.4226e+08 1.1515e+09 1018.98
## - PC29 1 8.0078e+08 1.2100e+09 1021.76
## - PC21 1 9.7120e+08 1.3804e+09 1029.14
## - PC13 1 1.1658e+09 1.5750e+09 1036.52
## - PC11 1 1.4845e+09 1.8937e+09 1046.84
## - PC5 1 2.2496e+09 2.6588e+09 1065.84
## - PC6 1 5.0428e+09 5.4520e+09 1106.06
## - PC1 1 3.4521e+10 3.4930e+10 1210.07
## - PC4 1 4.4955e+10 4.5364e+10 1224.71
## - PC3 1 9.4295e+10 9.4704e+10 1265.93
##
## Step: AIC=962.1
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 +
## PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 +
## PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + PC29 +
## PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 + PC39
##
## Df Sum of Sq RSS AIC
## - PC7 1 1.0512e+07 4.2751e+08 961.50
## <none> 4.1699e+08 962.10
## - PC19 1 1.6561e+07 4.3355e+08 962.28
## - PC39 1 1.8087e+07 4.3508e+08 962.48
## - PC32 1 2.3269e+07 4.4026e+08 963.14
## - PC28 1 2.3616e+07 4.4061e+08 963.19
## - PC34 1 2.5199e+07 4.4219e+08 963.39
## - PC35 1 2.5934e+07 4.4293e+08 963.48
## - PC10 1 2.7484e+07 4.4448e+08 963.68
## - PC8 1 2.7565e+07 4.4456e+08 963.69
## - PC16 1 5.0163e+07 4.6716e+08 966.46
## - PC25 1 5.6999e+07 4.7399e+08 967.28
## - PC20 1 5.9001e+07 4.7599e+08 967.51
## - PC30 1 6.6964e+07 4.8396e+08 968.44
## - PC38 1 8.1959e+07 4.9895e+08 970.15
## - PC37 1 9.4668e+07 5.1166e+08 971.56
## - PC17 1 9.6981e+07 5.1397e+08 971.81
## - PC18 1 1.1650e+08 5.3349e+08 973.90
## - PC33 1 1.8418e+08 6.0117e+08 980.59
## - PC26 1 2.0766e+08 6.2465e+08 982.73
## - PC22 1 3.0431e+08 7.2130e+08 990.79
## - PC9 1 3.2490e+08 7.4190e+08 992.36
## - PC31 1 3.7020e+08 7.8719e+08 995.68
## - PC12 1 3.8197e+08 7.9896e+08 996.51
## - PC15 1 4.7324e+08 8.9023e+08 1002.57
## - PC23 1 5.5779e+08 9.7479e+08 1007.65
## - PC24 1 5.6938e+08 9.8637e+08 1008.31
## - PC14 1 6.3476e+08 1.0517e+09 1011.91
## - PC2 1 7.4226e+08 1.1593e+09 1017.36
## - PC29 1 8.0078e+08 1.2178e+09 1020.12
## - PC21 1 9.7120e+08 1.3882e+09 1027.45
## - PC13 1 1.1658e+09 1.5828e+09 1034.80
## - PC11 1 1.4845e+09 1.9015e+09 1045.07
## - PC5 1 2.2496e+09 2.6666e+09 1064.01
## - PC6 1 5.0428e+09 5.4597e+09 1104.14
## - PC1 1 3.4521e+10 3.4938e+10 1208.08
## - PC4 1 4.4955e+10 4.5372e+10 1222.72
## - PC3 1 9.4295e+10 9.4712e+10 1263.93
##
## Step: AIC=961.5
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC8 + PC9 + PC10 + PC11 +
## PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 + PC20 +
## PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + PC29 + PC30 +
## PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 + PC39
##
## Df Sum of Sq RSS AIC
## <none> 4.2751e+08 961.50
## - PC19 1 1.6561e+07 4.4407e+08 961.62
## - PC39 1 1.8087e+07 4.4559e+08 961.82
## - PC32 1 2.3269e+07 4.5077e+08 962.46
## - PC28 1 2.3616e+07 4.5112e+08 962.51
## - PC34 1 2.5199e+07 4.5270e+08 962.70
## - PC35 1 2.5934e+07 4.5344e+08 962.79
## - PC10 1 2.7484e+07 4.5499e+08 962.98
## - PC8 1 2.7565e+07 4.5507e+08 962.99
## - PC16 1 5.0163e+07 4.7767e+08 965.71
## - PC25 1 5.6999e+07 4.8450e+08 966.50
## - PC20 1 5.9001e+07 4.8651e+08 966.73
## - PC30 1 6.6964e+07 4.9447e+08 967.64
## - PC38 1 8.1959e+07 5.0946e+08 969.32
## - PC37 1 9.4668e+07 5.2217e+08 970.70
## - PC17 1 9.6981e+07 5.2449e+08 970.94
## - PC18 1 1.1650e+08 5.4400e+08 972.99
## - PC33 1 1.8418e+08 6.1168e+08 979.56
## - PC26 1 2.0766e+08 6.3516e+08 981.67
## - PC22 1 3.0431e+08 7.3181e+08 989.60
## - PC9 1 3.2490e+08 7.5241e+08 991.15
## - PC31 1 3.7020e+08 7.9771e+08 994.43
## - PC12 1 3.8197e+08 8.0948e+08 995.25
## - PC15 1 4.7324e+08 9.0074e+08 1001.23
## - PC23 1 5.5779e+08 9.8530e+08 1006.25
## - PC24 1 5.6938e+08 9.9688e+08 1006.91
## - PC14 1 6.3476e+08 1.0623e+09 1010.47
## - PC2 1 7.4226e+08 1.1698e+09 1015.86
## - PC29 1 8.0078e+08 1.2283e+09 1018.60
## - PC21 1 9.7120e+08 1.3987e+09 1025.87
## - PC13 1 1.1658e+09 1.5933e+09 1033.17
## - PC11 1 1.4845e+09 1.9120e+09 1043.38
## - PC5 1 2.2496e+09 2.6771e+09 1062.23
## - PC6 1 5.0428e+09 5.4703e+09 1102.25
## - PC1 1 3.4521e+10 3.4948e+10 1206.10
## - PC4 1 4.4955e+10 4.5383e+10 1220.73
## - PC3 1 9.4295e+10 9.4722e+10 1261.94
summary(mod_ajustado)
##
## Call:
## lm(formula = Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC8 + PC9 +
## PC10 + PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 +
## PC19 + PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 +
## PC29 + PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 +
## PC39, data = df_pca)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6119.0 -1972.8 -442.5 1677.2 6079.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 391149.2 633.9 617.082 < 2e-16 ***
## PC1 7948.6 202.9 39.169 < 2e-16 ***
## PC2 1689.9 294.2 5.744 1.55e-05 ***
## PC3 -20651.3 319.0 -64.737 < 2e-16 ***
## PC4 -19536.4 437.1 -44.699 < 2e-16 ***
## PC5 -5201.2 520.2 -9.999 5.27e-09 ***
## PC6 8640.3 577.2 14.971 5.70e-12 ***
## PC8 733.6 662.8 1.107 0.282176
## PC9 -2792.0 734.7 -3.800 0.001210 **
## PC10 -889.0 804.4 -1.105 0.282864
## PC11 -7610.6 936.9 -8.123 1.34e-07 ***
## PC12 -4530.4 1099.5 -4.120 0.000582 ***
## PC13 8178.5 1136.2 7.198 7.75e-07 ***
## PC14 -6912.5 1301.4 -5.311 3.98e-05 ***
## PC15 6402.1 1396.0 4.586 0.000202 ***
## PC16 2178.5 1459.0 1.493 0.151822
## PC17 -3261.5 1571.0 -2.076 0.051699 .
## PC18 -3827.9 1682.3 -2.275 0.034651 *
## PC19 1670.9 1947.6 0.858 0.401630
## PC20 -3589.7 2216.8 -1.619 0.121857
## PC21 15742.9 2396.2 6.570 2.73e-06 ***
## PC22 9454.3 2570.8 3.678 0.001599 **
## PC23 15197.6 3052.3 4.979 8.33e-05 ***
## PC24 16195.3 3219.5 5.030 7.43e-05 ***
## PC25 6123.6 3847.4 1.592 0.127971
## PC26 13572.9 4467.8 3.038 0.006768 **
## PC28 -6166.5 6019.0 -1.025 0.318466
## PC29 41051.8 6881.3 5.966 9.66e-06 ***
## PC30 12587.3 7296.4 1.725 0.100730
## PC31 37646.4 9281.1 4.056 0.000674 ***
## PC32 -9955.8 9789.9 -1.017 0.321953
## PC33 -32279.3 11282.4 -2.861 0.009998 **
## PC34 15625.8 14765.4 1.058 0.303201
## PC35 20042.1 18668.1 1.074 0.296447
## PC37 102290.8 49868.9 2.051 0.054295 .
## PC38 -108663.4 56934.8 -1.909 0.071543 .
## PC39 -78836.9 87931.4 -0.897 0.381166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4743 on 19 degrees of freedom
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.9935
## F-statistic: 235.9 on 36 and 19 DF, p-value: < 2.2e-16
Encontramos um modelo com apenas \(36\) variáveis, grande parte delas é significante e apenas uma eliminação foi realizada. O coeficiente \(R^2\) é altíssimo, na ordem de \(99%\), indicando, a princípio, que uma grande parcela da variação quadrática total pode, de fato, ser explicada pelo modelo.
Porém, para termos uma noção geral dos resultados, precisamos olhar mais detalhadamente os resíduos.
Observemos os resíduos do modelo ajustado:
par(mfrow=c(2, 2))
mod_ajustado %>% aov %>% plot
anares <- resid(mod_ajustado)
ad.test(anares)
##
## Anderson-Darling normality test
##
## data: anares
## A = 0.38975, p-value = 0.3719
shapiro.test(anares)
##
## Shapiro-Wilk normality test
##
## data: anares
## W = 0.98017, p-value = 0.4828
Os p-valores são muito maiores que \(0.05\) (tanto para o teste de Anderson-Darling quanto para o teste de Shapiro-Wilk) e, assim, podemos adotar a hipótese de normalidade. Testando a homocedasticidade por meio do teste de Breusch-Pagan
bptest(mod_ajustado)
##
## studentized Breusch-Pagan test
##
## data: mod_ajustado
## BP = 34.627, df = 36, p-value = 0.5339
O elevado P-Valor não nos permite rejeitar a hipótese de homocedasticidade. Finalmente, testemos a autocorrelação por meio de um teste de Durblin-Watson:
dwtest(mod_ajustado)
##
## Durbin-Watson test
##
## data: mod_ajustado
## DW = 2.2876, p-value = 0.1314
## alternative hypothesis: true autocorrelation is greater than 0
Como o p-valor é maior que \(5 \%\). Assim, temos todas as hipóteses necessárias para uma boa regressão linear (já sabemos que não há multicolinearidade pois realizamos decomposição por componentes principais na transformação dos dados).
Nosso modelo cumpre os requisitos de uma boa regressão linear e podemos, agora, interpretar os resultados. Porém, para isso, precisamos verificar a influência de cada um dos clusters sobre cada um dos componentes principais. Isso será realizado na próxima seção.
Para interpretarmos o modelo, iremos, primeiramente, relembrar a sequência de passos adotada até aqui:
Assim, cada componente principal será explicado, em diferentes graus e percentuais, por diferentes clusters:
pr_comp_contrib <- pr_comp$rotation ** 2 %>%
apply(MARGIN=1, FUN=function(X) (X / sum(X))) %>% t %>% as.data.frame
pr_comp_contrib <- rownames_to_column(pr_comp_contrib, var='Cluster.Name')
Iremos demonstrar as contribuições para os valores principais mais significativos do modelo final:
summary(mod_ajustado)
##
## Call:
## lm(formula = Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC8 + PC9 +
## PC10 + PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 +
## PC19 + PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 +
## PC29 + PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 +
## PC39, data = df_pca)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6119.0 -1972.8 -442.5 1677.2 6079.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 391149.2 633.9 617.082 < 2e-16 ***
## PC1 7948.6 202.9 39.169 < 2e-16 ***
## PC2 1689.9 294.2 5.744 1.55e-05 ***
## PC3 -20651.3 319.0 -64.737 < 2e-16 ***
## PC4 -19536.4 437.1 -44.699 < 2e-16 ***
## PC5 -5201.2 520.2 -9.999 5.27e-09 ***
## PC6 8640.3 577.2 14.971 5.70e-12 ***
## PC8 733.6 662.8 1.107 0.282176
## PC9 -2792.0 734.7 -3.800 0.001210 **
## PC10 -889.0 804.4 -1.105 0.282864
## PC11 -7610.6 936.9 -8.123 1.34e-07 ***
## PC12 -4530.4 1099.5 -4.120 0.000582 ***
## PC13 8178.5 1136.2 7.198 7.75e-07 ***
## PC14 -6912.5 1301.4 -5.311 3.98e-05 ***
## PC15 6402.1 1396.0 4.586 0.000202 ***
## PC16 2178.5 1459.0 1.493 0.151822
## PC17 -3261.5 1571.0 -2.076 0.051699 .
## PC18 -3827.9 1682.3 -2.275 0.034651 *
## PC19 1670.9 1947.6 0.858 0.401630
## PC20 -3589.7 2216.8 -1.619 0.121857
## PC21 15742.9 2396.2 6.570 2.73e-06 ***
## PC22 9454.3 2570.8 3.678 0.001599 **
## PC23 15197.6 3052.3 4.979 8.33e-05 ***
## PC24 16195.3 3219.5 5.030 7.43e-05 ***
## PC25 6123.6 3847.4 1.592 0.127971
## PC26 13572.9 4467.8 3.038 0.006768 **
## PC28 -6166.5 6019.0 -1.025 0.318466
## PC29 41051.8 6881.3 5.966 9.66e-06 ***
## PC30 12587.3 7296.4 1.725 0.100730
## PC31 37646.4 9281.1 4.056 0.000674 ***
## PC32 -9955.8 9789.9 -1.017 0.321953
## PC33 -32279.3 11282.4 -2.861 0.009998 **
## PC34 15625.8 14765.4 1.058 0.303201
## PC35 20042.1 18668.1 1.074 0.296447
## PC37 102290.8 49868.9 2.051 0.054295 .
## PC38 -108663.4 56934.8 -1.909 0.071543 .
## PC39 -78836.9 87931.4 -0.897 0.381166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4743 on 19 degrees of freedom
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.9935
## F-statistic: 235.9 on 36 and 19 DF, p-value: < 2.2e-16
São as componentes: PC1, PC3 e PC4, que possuem significância abaixo de \(0.1\%\) no teste de Student. Plotando as contribuições:
pr_comp_contrib_p <- pr_comp_contrib %>%
gather(key='Component', value='Contrib', -Cluster.Name) %>%
arrange(desc(Contrib))
plot.Radar <- function(component_name) {
options(warn=-1)
df_to_plot <- pr_comp_contrib_p %>% filter(Component %in% c(component_name))
max_prob <- max(df_to_plot$Contrib)
size_x <- length(df_to_plot$Cluster.Name %>% unique)
seq_angle <- 360/(2*pi) * rev(pi/2 + seq(pi/size_x,
2*pi - pi/14,
len = size_x))
(ggplot(data=df_to_plot, aes(x=reorder(Cluster.Name, -Contrib),
y=Contrib, group=Component,
fill=Contrib)) +
geom_point(size=1) +
geom_bar(stat='identity') +
ggtitle('Análise de Componentes Principais' %>% paste(component_name, sep=': ')) +
geom_hline(aes(yintercept=0), lwd=1, lty=2) +
scale_y_continuous(limits=c(0, 1.35*max_prob)) +
coord_polar() +
theme_minimal() +
theme(axis.ticks =element_blank(),
axis.text.y =element_blank(),
axis.title=element_blank(),
axis.text.x=element_text(size = 6, angle=seq_angle))) %>%
return
}
plot.Radar('PC1')
plot.Radar('PC3')
plot.Radar('PC4')
Interpretando-se as componentes principais, obtemos:
Podemos então notar que o comércio internacional possui grande peso na determinação da emissão de CO2 pela França e, em expecial, o comércio de bens importados da América Latina.
Além disso, observamos uma participação relevante da produção de energia e da indústria bélica em tal variável. Logicamente, todos esses fatores estão correlacionados e representam faces diferentes de tal problema ambiental.
A análise das componentes principais na seção anterior nos forneceu alguns insights mas ainda não pareceu ser o suficiente para interpretarmos todas as variáveis explicativas (que são muitas) simultaneamente.
Por isso, propomos um método alternativo para mapear as variáveis relevantes do modelo: iremos gerar um wordcloud com as palavras com maior peso na determinação da saída (emissão de CO2 em KT).
Iremos criar um dataframe relacionando cada palavra contida na descrição das variáveis ao peso de tal variável na determinação da variável explicada. E como determinaremos tal peso?
Iremos utilizar a regra da cadeia para calcular a sensibilidade de cada variável na saída.
\(X_k = \sum_{i=0}^{i=N_{Clusters}} P_i.C_i\) e
\(C_k = \frac{\sum_{i=0}^{i=N_{NV_k}} V_{k,i}}{NV_k}\)
Onde:
Nesse caso temos:
\(\frac{\partial Y}{\partial V_{k, i}} = \frac{A_k \times P_i}{NV_k}\)
E esse será o peso de cada uma das variáveis, que será somado ao peso de cada uma das palavras contidas dentro da descrição da variável.
Com esses dados, iremos dispor os gráficos em uma nuvem de palavras (iremos utilizar apenas variáveis com p-valor abaixo de \(5 \%\) nessa análise).
concat_descript <- paste(df_in_meta$`Indicator Name`, collapse=' ') %>% toupper()
concat_descript <- gsub('[^[:alnum:] ]', '', concat_descript)
words_vec <- strsplit(concat_descript, split=' ') %>% unlist() %>% unique()
words_weight_vec <- zeros(length(words_vec)) %>% as.vector()
df_wordcloud <- data.frame(Word=words_vec,
Weight=words_weight_vec,
stringsAsFactors=F)
pc_names <- summary(mod_ajustado)$coefficients %>% rownames()
df_coef_sig <- summary(mod_ajustado)$coefficients %>% as.data.frame()
df_coef_sig$PC <- pc_names
df_coef_sig <- df_coef_sig %>% filter(`Pr(>|t|)` <= 0.05)
df_in_meta_grupos <- df_in_clusters %>%
inner_join(df_grupos, by=c('Cluster.Index'='Cluster.Index')) %>%
inner_join(df_in_meta, by=c('Indicator.Code'='Indicator Code'))
for (var_X in df_coef_sig$PC[-1]) {
var_A <- mod_ajustado$coefficients[[var_X]]
for (var_V in pr_comp$rotation[, var_X] %>% names()) {
weight_V <- pr_comp$rotation[[var_V, var_X]] * var_A
V_description <- (df_in_meta_grupos %>%
filter(Cluster.Cod==var_V) %>%
as.data.frame())$`Indicator Name.x`[[1]]
V_description <- gsub('[^[:alnum:] ]', '', V_description)
V_description_words <- V_description %>% toupper() %>% strsplit(split = ' ')
for (word_from_V in V_description_words %>% unlist) {
curr_val <- (df_wordcloud %>% filter(Word == word_from_V))$Weight[[1]]
setDT(df_wordcloud)[Word == word_from_V, Weight := curr_val + weight_V]
}
}
}
ignore_list <- c('', 'OF', 'PER', 'ANNUAL', 'NET',
'AND', 'OF', 'LCU', 'FROM', 'CONSTANT',
'TERMS', 'ADJUSTMENT', 'ON', 'TO', 'KM', 'KWH',
'SQ', 'DEC', 'IN', 'TOTAL', 'FINAL', 'BY')
df_wordcloud <- df_wordcloud %>%
filter(Weight != 0 & ! Word %in% ignore_list) %>%
mutate(Weight = abs(Weight)) %>%
mutate(Weight=scales::rescale(Weight, to=c(0, 20))) %>%
arrange(desc(Weight))
df_wordcloud
wordcloud(df_wordcloud$Word,
df_wordcloud$Weight,
colors = brewer.pal(8, 'Dark2'))
Óleo, produção, gás e combustíveis parecem exercer um primeiro plano nas palavras-chave que possuem mais peso na explicação das emissões de CO2.
População, bens e trocas também exercem um papel intermediário. Com relação a trocas, temos que, na seção anterior, encontramos relações fortes com diversas variáveis ligadas ao comércio internacional (importações e exportações).
Assim, nosso modelo, ao buscar as variáveis mais relevantes, encontrou três pilares determinantes na explicação das emissões de CO2:
E podemos concluir que a emissão de CO2 é um problema desafiante pois se encontra altamente correlacionada com termos-chave relacionados a desenvolvimento e crescimento econômico.
É por isso que a produção precisa ser contida e equilibrada com uma boa diretriz ambiental no conceito de desenvolvimento sustentável.
Para finalizar, iremos prever a emissão de CO2 para o ano de \(2014\) (último ano da análise) conforme o proposto no roteiro deste projeto. Para isso, iremos simplesmente utilizar a função predict.
Iremos também estimar as emissões para os anos de \(2012\) e \(2013\) pois também temos dados perdidos para a emissão de CO2 nesses anos:
df_predict <- predict.lm(object=mod_ajustado,
newdata=df_pca[c('2012', '2013', '2014'),],
interval='prediction') %>% as.data.frame()
df_predict
Plotando tais resultados:
df_pca$Year <-rownames(df_pca)
df_predict$Year <- c(2012, 2013, 2014)
df_pca_non_missing <- df_pca %>% mutate(Year = as.numeric(Year))
df_pca_non_missing <- df_pca_non_missing %>% filter(Year <= 2012 & Year > 2005)
ggplot() +
geom_line(data=df_pca_non_missing,
aes(x=Year, y=Y, group=1),
color='blue', stat='identity', group=1) +
geom_line(data=df_predict, aes(x=Year, y=fit),
colour='red',
linetype='dashed') +
geom_linerange(data=df_predict, aes(x=Year, y=fit, ymin=lwr, ymax=upr),
color='red', group=1) +
geom_point(data=df_pca_non_missing, aes(x=Year, y=Y),
shape=21,
colour='blue',
fill='white',
size=3,
stroke=1) +
geom_point(data=df_predict, aes(x=Year, y=fit),
shape=21,
colour='red',
fill='white',
size=3,
stroke=1) +
geom_point(data=df_predict, aes(x=Year, y=lwr),
shape=45,
colour='red',
size=15) +
geom_point(data=df_predict, aes(x=Year, y=upr),
shape=45,
colour='red',
size=15) +
labs(
title = 'Previsão para 2012, 2013 e 2014',
subtitle = 'Dados históricos em azul, previsões em vermelho',
x = 'Ano',
y = 'CO2 (KT) com Erro a 5% (Para previsões)'
) + theme_minimal()
O valor esperado para a emissão em KT no ano de \(2014\) é de \(339389.8\) em um intervalo compreendido entre \(327701.9\) e \(351077.6\) e a uma significância de \(5\%\).
Podemos ainda conferir tal resultado com o valor que, de fato ocorreu, para a emissão de CO2 na França nesse ano. Para isso, conferimos esse site: https://www.worldometers.info/co2-emissions/france-co2-emissions/, no qual encontramos uma emissão da ordem de \(320.703.500,00\).
O valor da emissão em \(2014\) se encontra muito acima daquilo que realmente ocorreu e a maior queda na emissão de gás carbônico na França ocorreu justamente nesse ano (\(\approx -9.0\%\)).
É plenamente possível que isso tenha ocorrido devido à implementação de políticas de redução na emissão de gás carbônico que não se encontram diretamente relacionadas às variáveis explicativas do modelo (porém, tal análise foge do escopo do presente projeto).
Em todo caso, isso demonstra que o modelo não é absoluto: os governos podem buscar uma crescente redução em tais emissões no intuito de implementar menos danos ambientais e isso pode ocorrer inclusive em variáveis não incluídas na base de dados aqui utilizada.
Os resultados são positivos para o caso do governo francês, observa-se uma tendêcia de queda.
Com base em cada uma das etapas, podemos tirar diversas conclusões: